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A  new  three-point  combined  compact  difference  (CCD)  scheme  is  developed  tor 
numerical  models.  The  major  features  of  the  CCD  scheme  are:  three  point,  implicit, 
sixth-order  accuracy,  and  inclusion  of  boundary  values.  Due  to  its  combination  ot 
the  tirst  and  second  derivatives,  the  CCD  scheme  becomes  more  compact  and  more 
accurate  than  normal  compact  difference  schemes.  The  efficient  twin-tridiugonal  (for 
calculating  derivatives)  and  triple-tridiagonal  (for  solving  partial  difference  equation 
with  the  CCD  scheme)  methods  are  also  presented.  Besides,  the  CCD  scheme  has 
sixth-order  accuracy  at  periodic  boundaries  and  fifth-order  accuracy  at  nonperiodic 
boundaries  The  possibility  of  extending  to  a  three-point  eighth-order  scheme  is  also 

included.  1)  I'WS  Ac*feniir  IV".. 


1.  INTRODUCTION 

The  grid  spacings  ( A.v,  Ay)  in  most  ocean  numerical  models  are  not  small.  For  example, 
a  global  ocean  model  is  considered  having  high  resolution  when  a  horizontal  grid  is  (1/8)  . 
approximately  14.5  km.  For  such  large  grid  spacing,  use  of  highly  accurate  difference 
scheme  becomes  urgent.  For  example.  McCalpin  [1|  used  fourth-order  differencing  to  re¬ 
duce  pressure  gradient  error  in  a -coordinate  ocean  models. 

The  trend  toward  highly  accurate  numerical  schemes  of  partial  differential  equations 
( PDE)  has  recently  led  to  a  renewed  interest  in  compact  difference  schemes.  Concurrently. 
Adam  [2],  Hirsh  [3],  and  Kreiss  [4]  have  proposed  Hermitian  compact  techniques  using 
less  nodes  (three  instead  of  five)  at  each  grid  point  to  solve  PDE.  Later  on.  as  pointed 
out  by  Adam  [5],  the  truncation  errors  are  usually  four  to  six  times  smaller  than  the  same 
order  noncompact  schemes.  Since  then,  much  work  has  been  done  in  developing  compact 
schemes  for  various  applications,  such  as:  an  implicit  compact  fourth-order  algorithm  |6|; 
a  fourth-order  compact  difference  scheme  for  nonunifomi  grids  [7];  fourth-order  and  sixth- 
order  compact  difference  schemes  for  the  staggered  grid  [S];  an  early  form  of  the  sixth-order 
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combined  compact  difference  scheme  |9[:  compact  finite  difference  schemes  with  a  range 
of  spatial  scales  (10):  and  an  ups'. inti  fifth-order  compact  scheme  [II].  These  schemes 
are  characterized  by  (a)  5-point  sixth-order,  (b)  much  lower  accuracy  at  nodes  adjacent  to 
boundaries,  and  let  no  requirement  on  PDE  to  be  satisfied  at  boundaries. 

Several  recent  work  emphasizes  on  the  improvement  of  boundary  accuracy.  For  hyper¬ 
bolic  system.  Carpenter er  al.  [12,  13]  introduced  a  simultaneous  approximation  term  (SAT > 
method  that  solves  a  linear  combination  of  the  boundary  conditions  and  the  hyperbolic  equa¬ 
tions  near  the  boundary.  This  method  provides  fourth-order  accuracy  at  both  interior  and 
boundary.  Under  the  assumption  that  the  derivative  operator  admits  a  summation-by  -parts 
formula  then  the  SAT  method  is  stable  in  the  classical  sense  and  is  also  time-stable.  For 
2D  vorticity  -stream  function  formulation.  E  and  Liu  [14.  15]  proposed  a  finite  difference 
scheme  with  fourth-order  accuracy  at  both  interior  and  boundary.  Question  arises:  can 
we  construct  a  scheme  (1 )  working  for  any  differential  equation  and  i2)  with  high-order 
accuracy  at  both  interior  and  boundary 

A  new  three-point  sixth-order  combined  compact  (CCD)  scheme  is  such  a  scheme  with 
the  following  features:  (a)  3-point  sixth-order,  (hi  comparable  accuracy  at  nodes  adjacent 
to  boundaries,  and  (c)  requirement  on  PDE  to  be  satisfied  at  boundaries.  Fourier  analysis  of 
errors  is  used  to  prove  the  CCD  scheme  as  having  better  resolution  characteristics  than  any 
current  ( uncompact  and  compact  l  scheme.  Two  impl  icit  solvers  for  the  CCD  scheme  are  also 
proposed  for  calculating  various  differences  (twin-tridiagonal  solver)  and  for  solving  PDEs 
(triple-tridiagonal  solver).  Furthermore,  we  use  the  one-dimensional  convection-diffusion 
equation  and  two-dimensional  Stommel  ocean  model  to  illustrate  the  application  of  the 
CCD  solvers  and  to  demonstrate  the  benefit  of  using  CCD  scheme. 


2.  CCD  SCHEME 
2.1.  Central  CCD  Algorithm 

Let  the  dependent  variable  f(x)  be  defined  on  the  interval,  0  <  .v  £  L.  Use  a  uniform 
grid.  0  =  .\  |  <  Xi  <  -tj  <  -  •  •  <  x\  <  xx  . i  =  L  with  a  spacing  h  —  .v,_i  —  .v,  =L/N.  Let 
the  dependent  variable  f(x  >  at  any  grid  point  x,  and  two  neighboring  points  .v,_  t  and  x,+t  be 
given  by  - 1  •  and  /._  i  and  let  its  derivatives  at  the  two  neighboring  points  .r,  _  i  and  , .  i 

be  given  by  //_, .  //I, . and  f'Y ,.  f"+i . //*J .  The  essence  of  the  CCD  scheme 

is  to  relate  f, .  /.’.  /" . f'k '  to  the  two  neighboring  points:  j  .  \ .  //l, .... ,  /,'* ’ 

;md  1.  y,_i .  /i+c  J:- 1- 


372 


CHU  AND  FAN 


and  to  compute  f," . //* 1  by  means  of  the  values  and  derivatives  at  the  two  neigh¬ 

boring  points.  Moving  from  the  one  boundary  to  the  other,  CCD  forms  a  global  algorithm  to 
compute  various  derivatives  at  all  grid  points.  In  this  paper  we  only  discuss  the  sixth-order 
CCD  scheme. 


2.2.  Local  Hermitian  Polynomial 

Let  //,  ( x)  be  a  local  Hermitian  polynomial  defined  on  the  closed  interval  v,+il, 
representing  the  variable  /  at  a,  and  /  and  its  derivatives  /'.  /"  at  the  two  neighboring 
points  .r,_j,  and  y.+|. 

)  =  /;-].  )  =  /(+|, 

,>  =  //_,.  //;<  v_,  )  =  | 

Expand  //,  (a  )  into  Taylor  series  in  the  neighborhood  of  with  sixth-order  accuracy 


H,  (a  )  =  H,  (Xj )  +  //  (  v,  ).v  -f 

Sf  ■  '  '•  * 


+ 


6! 


(2.3) 


The  seven  coefficients  in  (2.3)  are  detcnnined  by  the  seven  equations  in  ( 2.2 1. 


/<+)  -  ji- 1)-  t  +//.  i)+ — i  -  |) 


">■'  =  TSl/' 

//.'(a,)  =  ^j(/- 1  -2/,  +  /,_i)  -  ,  -  /,’.|i+  i +//!,)■ 

15 

4^ 

H 


15  15  ,  ,  3  „ 


-  -^Tt  -  2/,  +  I  -  —  I  -  //_,)-  j-j  l/'.  .  -  (.  j) 

"/V>  =  0(/;-t +  /;_,)  +  ^(/r;, 


(2.4) 


2/j' 

360 


(2.5) 


=  ^(/t-i  -2/;  +  /,.,>  -  +  /;:,). 

The  ill)  derivative  at  the  grid  point  .v,  is  approximately  given  by 

^  tf/'u,). 

Substitution  of  (2.5)  into  (2.4)  leads  to 

Z-i>  fi  ,6(^+i  A-i>  ~  tS  2/,(^+1  7781760^ 

^(Zi,-/;-,)-^/;:t+./;''i)+/1"  =  3^,-2/+/,  1>-^/;v 

(2.6) 
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which  ;trc  the  schemes  for  computing  the  first-order  and  second-order  derivatives  at  the  grid 
point  x, .  respectively.  Thus,  the  CCD  scheme  with  sixth-order  accuracy  can  be  written  by 


+(*i\  \  .  m  !l((6H\  _  fsiL\  \ 

16  \  \Sx)l+l  \ \SxJi  16  \  \Zx2  ),_J 


15  •  r 

—  ,  .,(/<•!  ft-  I  1 

16/t 

which  is  for  the  first  derivative  calculation,  and 


(2.7) 


=  I  ~  ~fi  +/(- 1) 


(2.8) 


which  is  for  the  second  derivative  calculation.  Comparing  (2.7 1  with  (2.1).  we  find  that  the 
parameters  in  (2.1)  for  the  sixth-order  scheme  should  be 


7  l 

Cf)  —  — .  p\  — - . 

16  16 


15  1  9 

to  = — .  cyi  = - .  fi- =  — .  =  3. 

X  "  X  4 


For  the  sixth-order  CCD  scheme,  the  truncation  errors  in  (2.6) 


1349 


-  ff'lf  s»  1.73*  I(r7'  h 


4 ,  f(  ‘ 1  t.b 


77X1760 


I 


20160 


/fV  as  4.9  *  10  V/S7/' 


are  quite  small. 

Another  benelit  of  using  ('CD  scheme  is  the  existence  of  a  global  Mcrmitian  polynomial 
w  ith  continuous  first-  and  second-order  derivatives  at  each  grid  point.  We  will  describe  it 
in  Appendix  1. 


2.3.  Error  Estimation 

We  compare  the  truncation  errors  between  the  CCD  scheme  with  current  generalized 
schemes  [10]  for  first-order  derivatives. 


f,  +  »{/.*,  *r  t",  ,  )  +  /l(./V+i+  =  « 


tl  T  I  / 1  ~  t 
2/i 


,Jit'  Ji-Z  ,/i-3  ,/i-J 
■  /> - — -  «‘ 


4/. 


6/r 


(2.9) 


and  the  second-order  derivatives. 

.  7)  —  I  -  -//  +  //  - 1  .fi-ri-lfi+fi-:  /,  +  .t  —  2/i  -r- 


/l- 


+.*- 


4/r 


+ 


9/t- 


(2.10) 


where  the  parameters  cz.  />. «,  />.  <  take  different  values  for  various  schemes  (Table  I ).  The 
comparison  of  truncation  errors  is  listed  in  the  last  column  in  Table  1.  We  find  that  the 
CCD  scheme  has  the  smallest  truncation  error  among  various  sixth-order  schemes.  For 
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TABLE  I 

Truncation  Errors  in  Various  Difference  Schemes  for  the  First 
and  Second  Derivative  Calculations 


Derivative 

approximation 

Parameter 

% 

Scheme 

a 

b 

c 

Truncation  error 

First 

,4, 

2nd-order  central 

0 

0 

i 

0 

0 

1 /<*/,-• 

(2.«j 

Standard  Padc  scheme 

1 

4 

0 

3 

3 

0 

0 

T' 

i2.ni 

6th-order  central 

0 

0 

3 

2 

-3 

~T 

1 

To 

36  * 

<2.f2) 

(«h -order  tri  diagonal 

1 

3 

0 

14 

T 

1 

9 

0 

4  x  2/  '/>" 

<2.  hi 

6th-order  pentadiagonal 

17 

57 

-) 

144. 

90 

57 

0 

n 

-too  1 

19  ‘  7!; 

(2.7) 

6th-order  CCD 

1 

/ 

/ 

/ 

/ 

1544  7! 

Second 

vO 

(2.13 1 

2nd-order  central 

0 

0 

1 

0 

0 

2  x  i/-/,-' 

i2.V5i 

Standard  Fade  scheme 

I 

10 

0 

6 

5 

0 

0 

“  x  V„- 
5  6!  ‘ 

[2.+^ 

(ith-ordcr  central 

0 

n 

3 

2 

—3 

5 

1 

ID 

72  x 

S'. 

(2.  rib 

hlh -order  (ridiagonal 

-> 

TT 

0 

12 

TT 

3 

TT 

0 

i  (  sr 

bill-order  pemadiagonul 

12 

97 

-1 

194 

120 

97 

0 

0 

—2672  1 

- -  x  —  f  ir 

(18) 

(ith-order  CCD 

/ 

1 

/ 

/ 

1 

_•)  x  2 

S'. 

example,  the  truncation  error  of  the  first  derivative  using  the  CCD  scheme  is  about  41.2 
times  smaller  than  using  the  sixth-order  central  scheme.  4.6  times  smaller  than  using  the 
sixth-order  tridiagonal  (compact  )  schemes,  and  6.0  times  smaller  than  using  the  sixth-order 
pentadiagonal  (compact)  scheme.  The  truncation  error  of  the  second  derivative  using  the 
CCD  scheme  is  about  36  times  smaller  than  using  the  sixth-order  central  scheme.  8.4  times 
smaller  than  using  the  sixth-order  tridiagonal  scheme  (compact),  and  13.8  times  smaller 
than  using  the  sixth-order  pentadiagonal  scheme  (compact).  Comparing  the  CCD  scheme 
with  the  second-order  central  difference  (SCD)  scheme  (most  commonly  used  in  ocean 
models),  truncation  errors  for  both  first  and  second  derivatives  are  more  than  four  orders  of 
magnitude  smaller. 

Another  good  feature  of  the  CCD  scheme  is  that  the  CCD  scheme  uses  the  same  formu¬ 
lation  at  all  grid  points  except  at  Lhc  boundaries,  where  some  additional  boundary  treatment 
is  formulated.  These  additional  schemes  at  the  boundaries  are  fifth-order  accurate  for  the 
PDE  w  ith  the  CCD  scheme  I  see  Section  5  ).  A  CCD  scheme  with  eighth-order  accuracy  w  ill 
be  presented  in  Appendix  2. 

3.  FOl  RIER  ANALYSIS  OF  ERRORS 

Fourier  analysis  of  errors  is  commonly  used  to  evaluate  various  difference  schemes, 
described  extensively  m  Swartz  and  Wcndroff  [16],  Oliger  and  Krciss  fl7|.  Vichnevctsky 


3-POINT  CCl?  SCHEME 


375 


and  Bowles  [18].  Roberts  and  Weiss  [  19].  Fromm  [20J.  Orszag  [21,  22|,  and  Leie  [10]. 
As  pointed  out  by  Lele  1 10].  Fourier  analysis  provides  an  effective  way  to  quantify  the 
resolution  characteristics  of  differencing  approximations. 

For  the  purpose  of  Fourier  analysis  the  dependent  variable  fix  \  is  assumed  to  be  periodic 
over  the  domain  [0.  L]  of  the  independent  variable,  i.e..  /i  =  /.v~i  and  h  =  L/ ,V .  The 
dependent  variable  may  decomposed  into  Fourier  series. 


*=.v,2 

/U)=  V 

it-=— iV/2 


(3.1 ) 


where  i  =  *f— 1.  It  is  convenient  to  introduce  a  scaled  wavenumber  »  =  2ir  khfL  — 
2nk/N .  and  a  scaled  coordinate  .v  =x/h.  The  Fourier  modes  in  tenns  of  these  are  simply 
expriu'.v  i.  The  exact  first-order  and  second-order  derivatives  of  (3.1 )  generate  a  function 
w  ith  exact  Fourier  coefficients 


••  i  w  » 

k  =  T/t. 


1  lowcver.  the  Fourier  coefficients  of  the  derivatives  obtained  from  the  differencing  scheme 
micht  not  be  the  same  as  the  exact  Fourier  coefficients,  i.e.. 


(&)/«/ A*  {ft  )ft>  =~ (  J  Jt- 


where  w‘  =  u  '(w )  and  it’"  =  u’"(u>)  are  the  modified  wavenumber  (both  real  numbers)  for 
the  lirst-order  and  second-order  differencing.  The  smaller  the  difference  between  the  exact 
and  modified  wavenumbers,  the  better  the  difference  scheme. 

According  to  Lelc  1 10|,  the  modified  wavenumbers  of  the  current  generalized  difference 
schemes  (2.9)  and  (2.10)  are 


w  (w)  = 


a  sin  w  -  't  sin  2 w  +  |  sin  3u’ 
I  -f  2a  cos  u:  +  Iff  cos  2  w 


(3.2) 


and 


i  r"(v>)=  y 


2a(l  -  cos u»)  +  ?[(1  -  cos2ir)4-  =^(1  -  cos3tr) 
1  2<y  cos  w  -|-2 cos  hr 


(3.3) 


respectively. 

For  the  CCD  schemes  (2.7)  and  (2.8).  the  modified  wavenumbers  a  '  and  w"  can  be 
calculated  joint  ly  as  follow  s: 


fu)  -  Y,he<,w" *“ 

t 

/'(.»,  =  Y'fc .  1 

L 

n.x)  =  y 


(3.4) 

(3.5) 

(3.6) 
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and 


k 

(3.7) 

i./r<v )  i ,,/  = 

k 

(3.8) 

fix  +h)  = 

k 

(3.9) 

fix -In  = 

■ 

(3.10) 

t 

(3.11) 

[/'(-«  =  V>T>/uc,r"',/MV 

k 

(3.12) 

lf\x  +  h)]r,  = 

X 

(3.13) 

1 

II 

M 

x 

(3.14) 

k 


Substitution  of  (3.4)— (3.14)  into  (2.7)-<2.8).  wc  have 


7,  ,,  ,  1  .  ,  15  . 

-[cos  u  +  1  ]«•  -t-  -  sin  v  (ir  —  sin  t/> 

8  8  8 


(3.151 


— (sin  w)w 
4 


I 

I - cos  if 

4 


(tf  )*  =  6{eos  w  —  I J. 


(3.16) 


Solving  (3.15H3.I6).  we  have 


ir'(ur)  = 

9  sin  if  [4  -t-  cos  if  | 

(3.17) 

24  -t-  20cos  w  -r  cos  2u 

11  "(Ml)  = 

!  8 1  -  48  cos  if  —  33  cos  2  if 
\  48  -)-  40  cos  w  -I- 2  cos  2u> 

(3.18) 

Among  various  difference  schemes,  the  modified  wavenumbers  of  the  tirst-order  differ¬ 
encing  w'  (Fig.  lai  and  of  the  second-order  differencing  u"  (Fig.  lb)  of  the  CCD  scheme 
are  closest  to  the  exact  wavenumber  w. 

In  multidimensional  problems  the  phase  error  of  first-order  differencing  scheme  appear 
in  the  form  of  anisotropy  [10.  18]. 


„  ,  ,  .(cos  ft)  mi  (uco.sft)  4-  (sinft)if  (irsinft) 

<(  Jfjlm.V)  =  w  (ir. «)/ir= - 

t  w 


(3.19) 


Figure  ic  shows  polar  plots  of  phase  speed  anisotropy  of  various  schemes  for  first  derivative 

approximations.  The  phase  speed  for  wavenumber  (magnitude)  w/n  =  jr,.  ^ . 5Tr  ii 

are  plotted.  Here,  we  also  see  that  the  CCD  scheme  shows  improvement. 
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(a) 

First  derivative  approximation 


1  1-5  ? 

Wavenumber 


(b) 

Second  derivative  apprcx-mation 


0  0.5  1  1.5  2  2.5  3 

Wavenumber 


(e) 


a  o  c  o  e  r 

Polar  plot  ot  phase  speed  anisotropy  tor  the  first  derivative  approximation; 
wavenumber; magnitude);  1/50,  5/50 . 45/50.  50/50 


FIG.  I.  Fourier  analysis  of  error  for  derivative  approximation;  tai  Second -order  centraJ  scheme;  lb)  standard 
Pade  scheme:  (c)  sixth-order  central  scheme;  id)  sixth-order  tridiacona)  scheme:  i.ei  sixth-order  pemadiaeonal 
scheme,  if)  combined  compact  scheme;  Igl  exact  differentiation. 


4.  CCD  FOR  DERI -VI  l\  E  CALCl  LATIONS 

The  previous  section  shows  that  the  sixth-order  3-pomt  CCD  scheme  is  more  accurate 
than  any  other  sixth-order  scheme  including  ordinary  compact  schemes.  Nevertheless,  since 
the  CCD  scheme  is  implicit  and  combines  computation  between  the  first-order  and  second- 
order  differences,  we  should  compute  /'  and  f"  jointly  and  globally. 

An  efficient  and  implicit  CCD  solver  is  designed  to  calculate  the  first-order  and  second- 
order  differences.  Since  CCD  is  3  3-point  scheme,  the  difference  calculation  at  v,  needs  to 
use  /.  /'.  and  f"  at  the  two  neighboring  points  v;  ,  and  \ ,  4 , .  At  the  two  boundaries  i ,  and 
.vy-ri.  some  specific  treatment  should  be  included  in  the  CCD  scheme. 


4.1.  Non-Periodic  Boundaries 

At  both  boundaries,  v  — i  |  and  v  =  .v\ .  i .  we  propose  a  fourth-order  one-sided  CCD 
scheme  instead  of  the  two-sided  scheme  to  keep  3-point  structure. 
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‘(8), '"‘(B).** 

ftr 

U*, 

J  -  +  />;/;  +  <>/:) 

(4.2) 

m  ...a)  .„( 

m 

1  =--(«(  /.v+l  +  *1  f\  +  <\f\  -1 ' 

iY  h 

(4.3) 

m 

\  =  rtds/v+i +<:/\ -i). 

V 

(4.4) 

\\  here 


«,=2.  /?,  =  -!.  «i|  =  —7/2.  b i=4.  fj  =  —1/2. 

C/2  =5.  $3  =  —6.  </;  =  0.  />;  =  —!  2.  £'3  =  3. 

Ai  the  boundaries,  ihe  frrsr-ordcr  difference,  represented  by  (4.1 1  and  (4.3).  has  a  truncation 
error  of  -|=  f"  h*.  The  second-order  difference,  represented  by  i4.2)  and  (4.4).  has  a 
truncation  error  of  —  fn'bi.  The  accuracy  at  both  boundaries  can  be  further  improved  to 
fifth  or  sixth  order. 

The  global  CCD  system,  consisting  of  (4.1 1  and  (4.2)  for  (  =  1,  (2.7)  and  (2.8)  for 

i  =  2, 3. 4 . .V.  and  (4.3)  and  (4.4)  for  i  =  .V  +•  1.  is  a  well-posed  system  since  it  has 

2(,V  4-  1  >  equations  with  2(N  +  1 1  unknowns:  (<)// &x ). .  f /6x:l,.  i  =  1. 2. 3 . .V. 

/V  +-  1.  We  may  write  the  2i.Y  1  >  equations  (4. 1 )— (4.4).  (2.7).  and  (2.8)  into  a  more 
general  form  (global  CCD  system). 


with 

«;<!)  =  /?,'(!)  =  =  f7v_,<3>  =  0.  y  =  1.2.  (4.6) 

representing  the  four  boundary  equations  (4. 1  m4.4  i.  Here,  j  =  I  corresponds  to  the  first- 
order  derivative  computation  (2.7).  and  j  —  2  corresponds  to  the  second-order  derivative 
computation  (2.8).  The  two  variables  .v /  and  s*  are  source  terms. 

Tlie  2( .V  +  1 1  x  2(iV  4-  1)  coefficient  matrix  of  (4.5 1  has  a  twin-tridiagonal  structure 
and  can  be  directly  solved  by  two  steps:  twin-forward  elimination  and  twin-backward 
substitution  (see  Appendix  3). 


4.2.  Periodic  Boundaries 

For  periodic  boundaries,  we  have 

tn  =  f,\ .  /i— /v-i-  7i.  =  /  \  •  f  i  ■  ff  =  h  —  f.v-c  14.7) 

Thus,  the  global  CCD  system,  consisting  of  (2.7)  and  (2.8)  for  /  =  1.2.3 . V.  is  well- 

posed  since  it  has  2.V  equations  with  2.V  unknowns:  (Sf/Sx), .  (<$-//<5.v;  i, .  /  =  1.2. 
3 . V.  The  coefficient  matrix  and  related  algorithm  are  listed  in  Appendix  4. 
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5.  CCD  FOR  SOLVING  FINITE  DIFFERENCE  EQUATIONS  (FDE) 

Any  PDE  discretized  by  the  CCD  scheme  (called  here  the  CCD  FDE)  can  only  be  solved 
globally  since  the  CCD  scheme  is  implicit.  Unlike  any  other  schemes,  the  CCD  FDE  solver 
requires  the  satisfaction  of  the  FDE  not  only  on  the  interior  points,  but  also  on  the  boundary 
nodes.  Benefits  of  such  a  treatment  are  to  decrease  the  truncation  errors  near  the  boundaries 
as  well  as  to  increase  the  global  accuracy.  Here,  we  propose  a  triple-tridiagonal  solver  for 
solving  CCD  FDE. 


5.1.  Nonperiodic  Boundaries 

Consider  a  one-dimensional  differential  equation. 

d  f  / 

tfi (a  )—  +a:i.t)—  -h  cZfd.v )/'(/)=  .v(.t ).  0  5  .r  <  L .  (5.1) 

dx  dx- 

with  general  boundary  conditions 

d\(x)f'(x)  -+-  J,,(.v)/(a  )  =  r(.v  l  at  v  =0:.v  =  /.,  (5.2) 

which  is  the  Dirichlet  boundary  condition  when  d„  =  1 .  </|  =0  and  the  Neumann  boundary 
condition  when  do  =  0.  d\  =  I . 

The  corresponding  FDE  can  be  written  as 

«un  (^)  +<»o (/>/  =.s, .  r  =  1.2 . V-l.  (5.5) 

and  the  boundary  conditions  become 

+  <(i)vT-,+4/W,=^  <5'4> 

Notice  that  we  applied  the  FDE  (5.3)  not  only  to  the  interior  points  but  also  to  the  two 
boundary  points  (  vi  and  a.v-m )•  At  each  interior  grid  n(xle  i  (2  <  /'  £  V)  we  have  three 
equations  [(5.3).  (2.7).  and  (2.8)]  with  three  unknown  variables  /,.  (tirf/Sx2),. 

However,  we  have  only  two  equations  [(5.3)  and  (5.4>|  at  both  boundaries  but  three 
unknowns:  /i.  (S//S.t)|.  (i2//'5.v’)j  for  the  left  boundary,  and  (8f/6x)\~\, 

|  for  the  right  boundary.  To  close  the  system  we  need  an  extra  condition 
for  both  the  left  and  right  boundaries. 

The  additional  boundary  conditions  are  obtained  by  constructing  a  new  fifth-order  poly¬ 
nomial. 


P (x )  =  Po  +  P\x  4-  Pi x1 4-  P}X}  4-  Pjx4  +  P<x\  (5.5) 

For  the  left  boundary,  the  six  coefficients  of  P(.r)  can  be  obtained  by 
/><.»,)  =  /,.  Pu2)=  f:.  Fir, )  =  /',.  P'(r,  )  =  /,'.  P,Ia2»  =  /2'.  = 

(5.6) 


The  additional  left  boundary  condition  with  fifth -order  accuracy  is  then  (Appendix  5' 


-32/2  +  /v)  =  0  (5.7) 


and  the  additional  right  boundary  condition  with  fifth-order  accuracy  is  written  as 


\ 


( 3 1  f\  *  i  —  32 /v  4-  fs  t )  =  0. 


15.8) 


Thus,  we  establish  three  equations  for  all  grid  points  i  interior  and  boundary  >  with  three 
unknowns  f] .  (8f/Sx)j.  ( S-f / 8.x2),.  i  =  1.2 . .V  +  |.  We  may  write  the  3<  V  -f  1 )  equa¬ 

tions  (2.7).  (2.8).  (5.3).  (5.4 1.  (5.7).  (5.8)  into  a  more  general  form  (global  CCD  FDE 
system). 


. 

+  />/(3)  ^ +c/(3)/,-i  —s'.  (5.9 


w  here  .•  =  1.2.3 . V  +  1  and  /  =  1 . 2.  3.  The  superscript  j  indicates  different  equations 

used  at  each  grid  point:  j  —  1  corresponds  to  FDE  (5.3).  j  =  2  corresponds  to  the  first-order 
derivative  calculation  (2.7).  and  j  =  3  corresponds  to  the  second-order  derivative  calculation 
(2.8).  For  all  the  interior  and  boundary  points,  the  coefficients  of  (5.9)  satisfy 


<1/(1)= 0/  (3)  =  A,-  ( I )  = b]  (3)  =  <‘  ( I )  =  t ,'  ( 3 ) = 0. 


(5.10) 


For  the  two  boundaries,  the  coefficients  of  (5.9)  satisfy 

o((l)=*f(l)=,c-f(l)  =  0, 

4+i<3)=H+i(3)=-4+,(3)=0.  ;  =  1.2.3.  (5.11) 

Thus,  the  coefficient  matrix  of  (5.9)  indicates  a  triple-tridiagonal  structure  and  can 
be  solved  in  two  steps:  triple-forward  elimination  and  triple-backward  substitution 
(Appendix  6). 


5.2.  Periodic  Boundaries 

17 

\  » 

For  periodic  boundaries  (4.9).  the  global  CCD  system  (5.9)  is  well-posed  since  it  has  3/V 

equations  with  3A  unknowns:  /,.  (8f/8x), .  (82f/8.x2),  .i—  1.2.3 . .V.  The  coefficient 

matrix  and  the  related  algorithm  are  listed  in  Appendix  7. 

6.  EXAMPLES 

The  CCD  scheme  proposed  here  is  a  three-point  scheme  w  ith  sixth-order  accuracy.  Usu¬ 
ally  a  three-point  scheme  te.g..  central  difference  scheme)  has  only  second-order  accuracy. 
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Two  examples  are  used  in  this  section  to  show  the  advantage  of  using  this  new  three- 
point  scheme.  Comparison  is  made  between  the  CCD  scheme  and  the  second-order  central 
difference  tSCDi  scheme  on:  (a)  truncation  error,  (h)  horizontal  resolution,  and  (cl  CPU 
time. 


6.1.  One-Dimensional  Convection-Diffusion  Equation 

Consider  a  one-dimensional  convection-diffusion  equation. 


,  i  irjr  , 

</(.v  iVv  +  b{.x ) — - ci.v)— —  —dix ). 

ax  dx~ 


0  <  x  •£  n. 


(6.1) 


with  the  boundary  conditions 


!//(())  =  0.  tft(.T>=0.  (6.2) 

If  the  coefficient  functions  in  (6.1 )  are  taken  as 

<*(.0  =  1.  h(x)~  I.  c(.v)  =  l,  dix)  —  cos.v  +  2 sin  v.  Os.'  <  .t.  (6.3) 


Eq.  (6.1 )  has  an  analytical  solution. 


if/un,(x)=  sin(x). 


(6.4) 


We  solved  (6. 1 1  numerically  with  both  CCD  and  SCD  schemes  under  various  horizontal 
resolutions,  and  we  recorded  the  CPU  time  (a  SUN  Sparc-20  was  used)  for  each  run. 
Comparing  the  numerical  results  with  the  analytic  solution  (6.4),  we  obtain  the  truncation 
errors  of  (he  two  schemes  for  the  given  resolution  (represented  by  number  of  cells).  We 
define  an  averaged  relative  error  (errav  >  by 


err„  — 


V  /^,,-^lA.tAy 

,  I'lVjA.vAy 


(6.5) 


Thus,  we  have  a  data  set  consisting  of  truncation  error.  CPU  time,  and  cell  number  for  the 
two  schemes. 

The  relationship  between  the  cell  number  UV)  and  err^  (Fig.  2a)  for  the  CCD  scheme 
(solid  curve)  and  the  SCD  scheme  (dashed  curve)  shows  that  for  the  same  err,v  the  cell 
number  would  be  much  smaller  in  the  CCD  scheme  than  in  the  SCD  scheme.  In  other 
words,  we  may  use  a  much  coarser  resolution  for  the  CCD  scheme  than  for  the  SCD 
scheme  if  the  same  accuracy  is  required.  For  example,  the  CCD  scheme  needs  only  18 
cells  when  crrJV  is  around  0.38  x  10  .  Flowever.  for  the  same  accuracy,  the  SCD  scheme 
requires  0400  cells  (sec  Tabic  2). 

The  relationship  between  the  CPU  time  and  the  averaged  relative  error  (Fig.  2b)  for  the 
CCD  scheme  (solid  curve)  and  the  SCD  scheme  (dashed  curve)  shows  that  for  the  same 
crrav  the  CPU  time  is  much  shorter  in  ihe  CCD  scheme  than  in  the  SCD  scheme. 

Such  striking  features  can  also  be  observed  in  Table  2.  When  the  relative  truncation  errors 
are  on  the  order  of  0.2  x  10  the  SCD  scheme  needs  .3600  grid  cells;  however,  the  CCD 
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10'8  10'7  1CT  10  *  10~* 

Average  error 


FIG.  2.  Comparison  between  the  CCD  anti  SCI)  schemes  in  one-dimensional  convection- -diffusion  equation: 
i  a)  cell  number  versus  average  error.  (biCPI  time  versus  average  error.  Here  solid  curves  denote  the  CCD  scheme 
am!  the  dashed  curves  represent  the  SC'D  scheme. 


scheme  requires  only  14  grid  cells.  The  CPU  time  is  also  more  than  an  order  of  magnitude 
smaller  using  the  CCD  scheme  (0.28  x  10  -1  s)  than  using  the  SCD  scheme  (0.32  x  10  1  m. 
The  ratio  of  CPU  between  using  SCD  and  CCD  schemes  (Ra).  called  the  CPU  ratio  here, 
is  around  24.2  when  ihe  truncation  errors  are  on  the  order  of  4.37  x  10  . 


6.2.  Stommel  Ocean  Model 

Stommel  [23)  designed  an  ocean  model  to  explain  the  westward  intensification  of  wind- 
driven  ocean  currents.  Consider  a  rectangular  ocean  with  the  origin  of  a  Cartesian  coordi¬ 
nate  system  at  the  southwest  comer  (Fig.  3).  The  x  and  v  axes  point  eastward  and  north¬ 
ward.  respectively.  The  boundaries  of  the  ocean  are  at  x  =0.  X  and  y  =  0.  b.  The  ocean  is 
considered  as  a  homogeneous  and  incompressible  layer  of  constant  depth  D  when  at  rest. 
Wien  currents  occur  as  in  the  real  ocean,  the  depth  differs  from  D  even where  b\  a  small 
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TABLE  2 


Comparison  between  the  CCD  and  SCD  Schemes  in  One-Dimensional 
Convection-Diffusion  Equation 


Error  range 

Features 

CCD 

SCD 

Ra 

0.36—0.83  x  1(1  - 

Cell  number 

7 

200 

Average  error 

0.3649  x  I0~* 

0.8292  x  10  4 

1.22 

CPI  time  (si 

0.0015 

O.OOI833 

0.27M>.35  x  II)  ’ 

Cell  number 

to 

1000 

Average  error 

0.2734  x  10- ' 

0.343  x  10  ' 

4.42 

C'Pl'  time  :s. 

0:002 

0.O08S33 

0-23—0.26  x  10' 

Cell  number 

14 

3600 

Average  error 

0.2395  x  10  » 

0,2577  x  I0-* 

11. 3 

CPC  time  (s) 

0.002833 

0.032 

0.37-0.38  X  10 

Cell  number 

18 

9400 

Average  error 

0.3747  x  10- 

0.3779  x  10  " 

24.2 

CPI  nmc  (s) 

0.0035 

0.08483 

perturbation.  Due  to  the  incompressibility,  a  streamfunction  ip  is  defined  by 


Hip  dip 

(/  =  -—,  t  =  — . 

dy  ox 

where  u  and  v  are  the  x  and  v  components  of  the  velocity  vector. 

Tlie  surface  wind  stress  is  taken  as  —Fco&ixy/h).  The  component  frictional  forces  are 
taken  as  -  Ru  and  —Rv.  w  here  R  is  the  frictional  coefficient.  The  Coriolis  parameter  f  is 
also  introduced.  In  general  it  is  a  function  of  y.  The  latitudinal  variation  of  f,  P  =  df/d\. 
is  called  the  /i -effect  in  the  ocean  dynamics.  Under  these  conditions  Stommel  derived  an 
equation  for  the  streamfunction  <// . 

(if1  3:  \  ,  34*  .  /,t  \ 

( — g  + — — =  — ysin  [  —  v  ).  (6.6) 

[ax2  ay- J  ax  {/>■  / 


yA 


b 


o 


* 


X 


FIG.  3.  Ocean  basin  dimensions  and  the  coordinate  system. 
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with  the  boundary  conditions 

(0.  >  )  =  4»(X.  y )  =  (jt .  0)  =  (x.  h)  =  0.  1 6.7 ) 

Here,  the  two  parameters  a  and  y  are  defined  by 

_  Dfl  _  F-r 
U~~R'  Y~~Rb' 

The  analytical  solution  of  (6.6)  vvith  the  boundary  conditions  (6.7)  is  given  by 

«P  =  —  —  j  sin  ^j-y^j(peu  +qeBx—  1),  (6.8) 

where 


a 

■-2  +  1 

a 2  /  rr  \  ~  cr 

i/t  +  UJ*  r=-2~\ 

'Hf)1 

(6.9) 

P  = 

(1 -<*)/(*■* -**■),  9= 

I  -  p. 

The  physical  parameters  are  -.elected  as |23] 

a  =  1 07  m.  h  =  2.t  x  1 0s  m.  D  =  200  m. 

F  =  0.3  x  10“7  nr  s“2.  R  =0.6  x  10  ’m s 

The  parameter  fi  is  taken  as  0  for  the  case  without  the  /(-effect  case,  and  it  is  taken  as 
10  11  m  1  s  1  for  the  case  w  ith  the  /(-effect  case. 


6.2.1.  Computational  Algorithm 

Use  a  uniform  grid.  0=.V|  <.v;<  •••  <.vy,  <  xN_|  =},.  and  0=V|  <  y-  <  •••  < 
yv  <.v.v,+i  —  h  with  grid  spacing  A.v  =  —  .v,  =  a /AT,  and  Ay=y,+i  -y, 

For  simplicity  and  no  loss  of  generality,  we  assume  that  the  cell  number  in  both  the  a  and  y 
directions  are  the  same.  /V,  =  N\.  =  ,V.  The  alternating  direction  implicit  (ADI)  method  is 
used  for  sols  ing  FDE.  The  iteration  k  to  A  —  I  can  he  separated  into  two  parts:  (a)  iteration 
along  the  r-axis  to  obtain  "intermediate  variables"  .  and  (S^V/Sx2)’ ,. 

1 V  9  f(s*Y  -  f—)*  ) 

s  V  iv2  A.,_,  +  8A.V  [\Sy  _  V  Sy  JLj.J 

—( f—Y  l  (-Y  \ .  _^l(  f^Y  ~  f—Y  \ 

16 /,+!.,  /  +  \&x),.j  16VV^-JyiAl/  V'*2/,-!.,/ 


;+ 1 


(6.10) 


(6.11) 


-l  (f—Y  ~ f—Y  ]--((— Y  (—Y  \ 

SA-t  ^  V  M  ) V  IS-  ,  I  8  v  \  S.\2  +  {  (i.v-  )  ,  / 

/  S2'P  \  ’  I 

-I-  — -  -i — (Apr ,  _  sxp-  +  vp-  )=o 

V  * v-  y , .  aa-  --1'  **v**M-r 


(6.12) 
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and  (b)  iteration  along  the  v-axis  to  obtain  variables  at  the  next  iteration  k  4-  I.  4'“  '  . 

(i4»/5.v)ji,,.and(32'|//3.T-)i;1. 


15  1 

X  2  Ay 


_  0/‘  +  ' 

i 


)=0 


(6.14) 
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(6.15) 


Such  an  iterative  process  stops  when  the  correction  at  the  iteration  k  - hi. 


corr 


(t+u 


V  vjP  • 

i  i 


k~  i 


4''  A  X  A  v 


V1 


-  . 


'K  ,  |  A.v  A  v 


(6.16) 


is  smaller  than  10 'f\ 


6.2.2.  Case  I:  Without  the  Effect 

The  condition  fi  =  0  leads  to  a  =  0  in  (6.6).  The  analytical  solution  of  (6.6)  becomes 


which  is  depicted  in  Fig.  4. 

We  solved  (6.6)  numerically  with  both  CCD  and  SC'D  schemes  under  various  horizontal 
resolutions,  and  we  recorded  the  CPU  time  (a  SUN  Sparc-20  was  used)  for  each  run. 
Comparing  the  numerical  results  with  the  analytic  solution  (6.17).  we  obtain  die  truncation 
errors  of  the  two  schemes  for  various  resolutions  (represented  by  the  number  of  cells). 

The  relationship  between  \  and  errttv  (Fig.  5a)  for  the  CCD  scheme  (solid  curve)  and 
the  SCD  scheme  idashed  curve)  shows  that  for  the  same  errav  the  cell  number  (A  i  would 
tie  much  smaller  for  the  CCD  scheme  than  for  the  SCD  scheme.  Tins  is  to  say  that  we  may 
use  a  much  coarser  resolution  for  the  CCD  scheme  than  for  the  SCD  scheme  for  the  same 
accuracy.  The  relationship  between  the  CPU  time  and  the  averaged  relative  error  (Fig.  5b) 
for  the  CCD  scheme  (solid  curve)  and  the  SCD  scheme  (dashed  curve i  shows  that  for  the 
same  err.,*  the  CPU  time  is  much  shorter  in  the  CCD  scheme  than  in  the  SCD  scheme. 

Table  3  lists  err;iv.  cell  number.  CPU  time  for  the  two  schemes,  and  CPU  ratio  i  fia ). 
When  the  relativ  e  truncation  errors  are  on  the  order  of  0.68  x  10  1 .  the  SCD  scheme  needs 


CPU  lime  (s)  Average  error 


%  10* 


FIG.  4.  Stieamfunciion  ( nr/M  obtained  from  Siommel  ocean  model  with  beta  —  0. 


FIG.  5.  Performance  of  the  CCD  and  SCD  schemes  in  Siommel  ocean  model  (beta  =01:  (a)  average  error 
versus  cell  numher  in  the  SCD  scheme:  (b)  average  error  versus  cell  number  in  the  CCD  scheme:  ic)  CPF  time 
versus  cell  numher  in  the  SCD  scheme:  icl  i  CPU  time  versus  cell  number  in  the  CCD  scheme. 
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TABLE 3 

Comparison  between  the  CCD  and  SCI)  Schemes  in  Slominel  Ocean  Model  (beta  =  0> 


Error  range 

Features 

CCD 

SCD 

Ra 

0.86-0.0  x  10  * 

Cell  number 

9x9 

50  x  50 

Average  error 

0.866  x  10  1 

0:894  x  10  * 

27.0 

CPI'  time  is) 

3.10 

83.8 

0,76-0.77  x  10  1 

Cell  number 

tO  x  10 

lot)  X  too 

Average  error 

0.766  x  10  * 

0.761  x  for* 

271.7 

CPC  time  (s) 

4.6 

1250 

0.68-0.69  x  10 

Cell  number 

14  x  14 

150  x  150 

Average  errni 

0.685  x  10"* 

0.68  x  10  * 

356.8 

CPC  time  ts) 

16-7 

5780 

22.500  grid  cells:  however,  the  CCD  scheme  requires  only  196  grid  cells.  The  CPU  ratio 
between  using  SCO  and  CCD  schemes  ( Ra)  is  356.8. 

6.2.3.  Case  2:  With  the  ft -Effect 

For  this  ease,  fl  =  10“ 1 1  rn_  !s  1  is  used.  The  analytical  strcamfunction.  is  plotted 
in  Fig.  6.  Wc  solved  (6.6)  numerically  with  both  CCD  and  SCD  schemes  under  various 
horizontal  resolutions,  and  we  recorded  the  CPU  time  (a  SUN  Sparc-20  was  used)  for 
each  run.  Comparing  the  numerical  results  with  the  analytic  solution  (6.8).  wc  obtain  the 
truncation  errors  of  the  two  schemes  for  various  given  resolutions  (represented  by  the 
number  of  cells). 

The  relationship  between  N  and  errav  (Fig.  7a)  for  the  CCD  scheme  (solid  curve)  and  the 
SCD  scheme  (dashed  curve)  shows  that  for  the  same  err*  the  cell  number  t.V )  would  be 
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FKf.  6,  Streamfunetion  (m7s>  obtained  from  Stommel  ocean  model  with  heta  =  10  :  m  < 
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Celi  number 


Fl«.  7.  Performance  of  the  CCD  and  SCD  schemes  in  Stommel  ocean  model  (beta  -  Id  m1  >  1  r 
■  a  )  average  error  versus  cell  number  inlhc  SCD  scheme.  <b:  average  error  versus  cell  number  in  tire  CCD  scheme. 
Id  CPC  time  versus  cell  number  in  the  SCD  scheme;  id.  CPI  time  versus  cell  number  in  Use  CCD  scheme. 


much  smaller  in  the  CCD  scheme  than  in  the  SCD  scheme.  The  relationship  between  the 
CPU  time  and  the  averaged  relative  error  (Fig.  7b >  for  the  CCD  scheme  (solid  curvet  and 
the  SCD  scheme  (dashed  curve)  shows  that  for  the  same  err,,  the  CPI  time  is  much  shorter 
in  the  CCD  scheme  than  in  the  SCD  scheme. 

Table  4  lists  err,v.  cell  number.  CPU  lime,  and  Ru  for  the  two  schemes.  When  the  relative 
truncation  errors  are  on  the  order  of  0.73  x  10  ’.  the  SCD  scheme  needs  22.500  grid  cells: 
however,  the  CCD  scheme  requires  only  729  grid  cells.  The  CPI  ratio  between  using  SCD 
and  CCD  schemes  t  Ra)  is  254.87. 


7.  CONCLUSIONS 

( 1 1  From  this  studs,  it  can  be  stated  that  the  three-point  sixth-order  CCD  scheme  is  a 
promising  highly  accurate  method  for  both  derivative  computation  and  FDE  solutions.  The 
advantage  of  this  scheme  is  the  existence  of  a  global  sixth-order  polynomial  which  not  only 
satisfies  the  FDE  at  all  the  grid  nodes  including  boundary  points  but  also  the  boundary 
conditions. 
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TABLE  4 

Comparison  between  the  CCD  and  SCO  Schemes  in  Stommcl 
Ocean  Model  (beta  =  10  11  m  1  s  1 1 


Error  range 

Features 

CCD 

SCD 

Rn 

0.2<H}.24  x  It) 

Cell  number 

14  <  14 

50  x  50 

Average  error 

0.23ft  x  in  - 

0  204  x  10- ; 

1.98 

CPI  timcM 

8.12 

16.1 

0.22-0.24  <10 

Cell  number 

l‘>  x  10 

150  x  150 

Average  error 

0.238  x  10  ' 

0.225  x  10 

78.79 

CPL  lime  tst 

14.9 

1174 

0.73-0.74  <io  • 

Cell  number 

27  x  27 

250  x  250 

Average  error 

0.73  x  HI-1 

0.735  k  1(1  ‘ 

25a.  87 

CPL'  tune  isi 

33.9 

SMO 

(2)  Fourier  analysis  shows  that  the  CCD  scheme  has  the  least  error  among  other  same 
order  schemes,  including  the  normal  compact  scheme.  Also,  the  CCD  scheme  has  the 
smallest  truncation  error  among  various  sixth-order  schemes.  The  truncation  error  of  the 
first  derivative  using  the  CCD  scheme  is  about  4 1 J  times  smaller  than  using  the  sixth-order 
central  scheme.  4.6  times  smaller  than  using  the  sixth-order  tridiagonal  (compact)  scheme, 
and  6.0  times  smaller  than  using  the  sixth-order  peniadiagona!  (compact)  scheme.  The 
truncation  errorof  the  second  derivative  using  the  CCD  scheme  is  about  36  times  smaller  rhan 
using  the  sixth-order  central  scheme.  X.4  times  smaller  than  using  the  sixth-order  tridiagonal 
scheme  (compact),  and  1 3.8  times  smaller  than  using  the  sixth-order  pentadiagonal  scheme 
(compact).  Comparing  the  CCD  scheme  w  ith  the  second-order  central  difference  (SCD.i 
scheme  (most  commonly  used  in  ocean  models),  the  truncation  errors  for  both  tirst  anil 
second  derivatives  are  more  than  four  orders  of  magnitude  smaller. 

(3 1  For  periodic  boundaries,  the  CCD  scheme  has  sixth-order  accuracy  ai  all  grid  points 
including  boundary  nodes.  For  nonperiodic  boundaries,  the  CCD  scheme  has  sixth-order 
accuracy  at  all  interior  grid  points,  fourth-order  accuracy  in  the  derivative  computation,  and 
fifth-order  accuracy  in  the  FDF  solutions  at  the  boundary  nudes. 

(4)  Both  twin-tridiagonal  and  triple-tridiagonal  techniques  are  proposed  for  the  CCD 
scheme  for  calculating  derivatives  and  solving  FDEs. 

(5)  Two  examples  (the  convection-diffusion  model  and  the  Stommcl  ocean  model)  show 
striking  results  (great  reduction  in  truncation  error  and  CPI  time),  which  may  lead  to  a 
wide  application  of  the  CCD  scheme  in  computational  geophy  sics. 

(6)  Future  studies  include  applying  the  CCD  scheme  to  nonuniform  and/or  staggered 
grid  systems,  as  well  as  designing  even  higher  order  schemes  such  as  an  eighth-order  CCD 
scheme. 


APPENDICES 

Appendix  l:  Global  Herinitian  Polynomial 

The  tirst-order  and  second-order  CCD  differences  are  obtained  implicitly  and  globally 
by  the  two  joint  equations  (2.7)  and  (2.N).  A  twin-tridiagonal  technique  was  developed  to 
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compute  /'  and  f"  at  all  grid  points.  As  soon  as  the  global  tirst  and  second  differences 
are  obtained,  the  higher  order  (k  —  3. 4.  5.  6)  differences  can  easily  be  calculated  locally 
with  (2.3). 

Since  the  CCD  scheme  is  solved  globally,  the  neighboring  local  Hermitian  polynomials 
should  satisfy 

//>,)  =  //;_,(.»,)  =  //;_ ,  (a,) 

Hit *.)  =  h;_x(x,)  =  H^  |(.c,) 

A  global  polynomial  //A. f.v)  can  be  defined  by 

H., (.V )  =  //:(.v).  a  — A  |  <  -V  <  A;. 

H.Ax)  =  +  (  I  —  w,  ty  <  x  <  .v,+i  (<  =2.3 . n—  1). 

//,(.v)  =  x„  <x<  .vn+i  =b, 

where  at,  (/  =  2. 3 . n  —  1 )  are  the  local  weighting  factors.  Notice  that  no  matter  what 

value  of  to,  is.  the  global  polynomial  H.Ax)  always  has  continuous  first-  and  second-order 
derivatives  at  the  point  v ... 

Hjx,)  =  Hjx,  -0)-Hjx,  +0) 

H”(x, )  =  H;\x,  -  0)  =  W;<a,  +  01. 

The  weighting  factors  are  recommended  to  be  0  <  «>,  <  I .  If  only  the  first-order  and  second- 
order  derivatives  are  computed,  we  may  use  at,  =  1  2  for  simplicity.  It  is  also  possible  to 
optimize  w,  by  minimizing  the  discontinuity  properties  of  the  high-order  <4  >  3)  deriva¬ 
tives  at  the  node  points.  As  soon  as  the  global  polynomial  HAx)  is  established,  we  can 
calculate  all  the  derivatives  and  integrate.  Since  the  values  of  at,  do  not  affect  the  first- 
order  and  second-order  derivatives,  we  will  not  discuss  here  the  effect  of  <w(.  This  pa¬ 
per  focuses  only  on  the  first-order  and  second-order  differentiation  of  the  second-order 
PDE. 

Furthermore,  a  higher  order  (higher  than  sixth-order)  three  points  CCD  scheme  can  also 
be  defined.  See  Appendix  2  for  description. 


Appendix  2:  Eighth-Order  CCD  Scheme 

The  eighth-order  CCD  scheme  relates  f.  f'.  f".  f  "  to  the  two  neighboring  points: 
/'-I-  //  i"  f,"- 1’  fi  I  and  /,+(,  and  solves  for  A  local 

Hermitian  poly  normal  H,ix )  is  defined  on  the  closed  interval  |  .v;  .  v, , ,  |  bv 

Hax)  =  HAXi)-^H,  toi)x  +  .v-r-^ — v  +  — , —  v  -p  — ,c 

H*'<Xi)  6  w;7,(.»y)  ,  Hf'lx,)  , 

H  6i  "*  7!  '  8!  *' 
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with 


// m.  i)  =  /"  - 1 ,  h, (.V; )— f.  ii  (a  w;«.v, «;cv, + 1 ) = /'. 


+1- 


r  9  .  t  —  f:l3> 


The  nine  parameters  are  determined  by 


H,(x.)  =  f, 


f  +  ~(.Ct  -I  -b/.-l) 

^(./i  +  l  -2/,  +/:-!>-  —  //_|)  + +  //-l) 

-4  (/)*-/,-!) 


:  (  v  i 


^  -/<-«)  +  +  //-I*- £(/?-!  "//-I’ 

+  (//ii +  //-t) 


w.,4>  = 


72  1 83  39 

—2/,  +/.-i)  +  ^]r(//T|  —  / - 1 >  —  +/"  i> 


//;s,(.v() 


^</+l  —/«-!>—  +/-!>+  jjjUfll  -/-!> 

^  '  f(3l  ,  /•(  3)  i 


w;6>(.v,) 


4/t 
1440 

~7i*~ 


l/+i  -2/,  +/,_i)- i  +  TT(//+i  +  /  - 1 ' 


//J 


45  )  ft?)  _  ,13), 

2pl/i+l  /i-tl 


W|»  3  -  - 


W  '(.v,)  = 


2575  .  ,  1575  ,,  ,  315 

</-.  - /i-.)r+  2jT</:+,  +/-«)  -  -srU+1  -  //I.) 


2/t7 

105 

lii1 


///S'(.v,» 


21060  „  .  ^  13860  „  3780 

'(/  fil  —  -/  +  h -I  >  :  ■>/)'  +  /-|) 


A* 


/i* 


592 


CHI.  AN  D  FAN 


The  k th  derivative  at  the  grid  point  .r,  i>  approximated  by 


/'*'(. v,)  =  ///* '  ( A", ).  *=1.2 . 8. 


Therefore,  the  first-order  derivative  at  grid  point  v,  is  computed  by 

+  fi-  i 1  +  A  ~  ^(/,it -/.-i>  + +//  'i) 

_  35  I  427  ,y(  h* 

~  \62fi"'  ~  1737  ^  8f 

the  second-order  derivative  at  grid  point  a,  is  computed  by 

29  .  ..  _.  5  ...  . .  h 


1  (ih 


48 


=  4^,-2/, +  f,-0  +  ±/r£: 

and  the  third-order  derivative  at  grid  point  a  is  computed  by 

-^«r.i+ tu  >  ♦  £<&,  T  cc,  +  ) + if 


_L  1 1.  _  f  i  -  f'yt  -_i_ 

8  2h-J"'  '  1  16212''  6! 


1357 w/»t’ 


Appendix  3:  Nonperiodic  CCD  Calculation 

Twin-forward  elimination/backward  substitution  scheme  is  designed  to  solve  global  CCD 
system  (4.5)  with  boundary  conditions  (4.6).  The  2(  X  —  I  >  x  2<  V  +  1 )  coefficient  matrix  of 
(4.5 1  has  a  tv.  in-tridiagonal  structure  and  can  be  directly  solved  by  two  steps:  twin-forward 
elimination  and  twin-backward  substitution. 

A.3.1.  Twin-Forward  Elimination 

The  twin-forward  technique  is  used  to  transform  the  tw  in-tridiagonal  coefficient  matrix 
into  a  twin-diagonal  coefficient  matrix  by  eliminating  the  four  parameters.  <  I ).  />’(!  >. 
<i;i  I )./»:( 1 )  at  each  grid  point  (Fig.  8).  At  the  left  boundary  u  =  1 ).  these  four  parameters 
are  already  absent. 

If  the  four  parameters  at  grid  node  i  are  eliminated,  it  is  easy  to  use  (4.5 1  to  eliminate 
aj+ 1  ( I ).  «r.  |  ( I ».  b)  .,(11.  />-. ,  ( 1 )  at  grid  point  /  4- 1 .  This  process  continues  until  reaching 
the  right  boundary.  The  coefficient  matrix  of  the  global  CCD  system  becomes  twin-diagonal. 
Figure  9  shows  the  structure  of  the  coefficient  matrix  after  twin-forward  elimination,  where 
the  shadowed  area  shows  the  eliminated  elements. 


A.3.2.  Twin-Backward  Substitution 

The  twin-backward  substitution  technique  is  used  to  obtain  both  ( Sf/Sx ),  and  (<5:  f/8x~), 
from  known  (Sf/Sx)  ,  i  and  (82f/8x:),+ 1.  After  the  twin-diagonal  coefficient  matrix  has 
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K1G.  8.  Siruciuiv  of  the  CCD  coefficient  malm  for  nonperiodic  boundaries. 


FKi.  9.  The  twin -forward  ciiminauon  of  the  CCD  cofftcicm  matrix  for  nonperiodic  boundaries.  Here  S 
denotes  eliminated  coefficients. 
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been  established,  the  global  CCD  system  (4.5)  becomes  two  equations  with  two  unknowns 
at  the  right  boundary  (Av-M ). 


Solving  this  set  of  two  algebraic  equations,  we  obtain  i  affix  )\  ..i  and 
The  substitution  procedure  starts  from  the  second  right  point  tvv).  The  first- and  second- 
order  differences  (8f/&x);  and  (S:f/ix2),  are  computed  front  substitution  |/  =  ;V, 
.V  —  1 . It: 


Appendix  4:  Periodic  CCD  Calculation 

TIte  structure  of  the  periodic  CCD  matrix  is  shown  in  Fig.  10.  Similar  to  nonperiodic 
boundaries,  we  construct  another  form  of  twin-forward  elimination  and  twin-backward 
substitution  procedures  for  periodic  boundaries.  Figure  1 1  show  s  the  structure  after  the  twin- 
forward  elimination  procedure,  where  the  shadowed  areas  mean  the  eliminated  elements. 

Appendix  5:  Fifth-Order  Accurate  Nonperiodic  Boundary  Conditions 

Consider  the  left  boundary  with  uniform  grid  A.v  =  It.  Let  vt  be  the  left  boundary  node: 
let  v;  and  v?  be  the  first  and  second  neighboring  nodes.  Expanding  the  dependent  variable 
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H<>.  1 1.  The  twin- 1  orvard  elimination  of  1  he  CCD  coefficient  matrix  for  periodic  boundriev  Here  E3  denotes 
eliminated  coefficients. 


h  f _ I  )- 

/ Ut)  -  fix 2J+  Vi- i-/a'(jr:)/|^0(/r) 

Jt=i 

**  I 

/(  a-3)  =  f(xi)  +  2I  n/^ta^  +  ocA7) 

*=i  *■ 

*  <_h* 

/*Ut )  =  /'(.t2)  +  ]T  /“  •  h(.v2)/;*  -  <?<//’) 

<-i  ■ 

/"(.V,)  =  53  m^/,4i2*l.*a)^+0(AS) 

4=1 


which  lead  to 

14/'(.v,)+  l6/'(.v:)  -2/''(.vt)/t  -  4/"Ui'Mr  +  }.(3l/(Jvr)-32/(x2).+  /(jr3» 

h 

=  l&J*'(X2)  +  Oih*). 


Therefore,  the  nonperiodic  boundary  condition 


has  fifth-order  accuracy. 
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Appendix  6:  Nonperiodic  CCD  FDE  Solution 

Tlie  triple-forward  elimmation/backsvard  .substitution  scheme  is  designed  to  solve  global 
CCD  FDB  system  (5.9)  with  boundary  conditions  (5.10).  The  3(  A  +  I)  x  3  (Ai  -  I)  coef¬ 
ficient  matrix  of  1 5.9)  has  a  triple-tridiagonal  structure  and  can  be  directly  solved  by  two 
sieps:  triple-forward  elimination  and  triple-backward  substitution. 


A.6. 1 .  Triple-Forward  Elimination 

The  triple-forward  technique  is  used  to  transform  the  triple-tridiagonal  coefficient  matrix 
into  a  triple-diagonal  coefficient  matrix  by  eliminating  the  six  parameters,  a)  ( I),  a~(  I ). 
h\  ( 1 1.  Ik(  I ).  <■'  1 1 ).  c2l  1 )  at  each  gnd  point  ( Fig.  12).  At  the  left  boundary  (/  =  1 ).  these  six 
parameters  are  already  absent. 

If  the  six  parameters  at  grid  node  i  arc  eliminated,  it  is  easy  to  use  (5.9)  to  eliminate 
i/.Ljfl  ).  I ).  b}+ ,( 1 ).  h;. , ( 1 ),  <•,'(  1 },  r;(  1 1  at  grid  point  / -f  I.  This  process  continues 
until  reaching  the  right  boundary.  The  coefficient  matrix  of  the  global  CCD  F'DE  system 
becomes  triple-diagonal.  Figure  1 3  shows  the  structure  of  the  coefficient  matrix  after  triple¬ 
forward  elimination,  where  the  shadowed  area  shows  the  eliminated  elements. 


A.6.2.  Triple-Rat  kward  Substitution 

The  triple-backward  substitution  technique  is  used  to  obtain  /, .  (Sf/Sx ), .  and  I f/Sx2), 
from  known  /1+] .  (Sf/&x >.  + 1 .  and  (<$:  f/S.x2  i,„i  .  After  the  triple-diagonal  coefficient  matrix 
has  been  established,  the  global  CCD  system  (5.9)  becomes  three  equations  with  three 
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FI< >.  12.  Structure  of  the  CCD  coefficient  mams  for  FDE  with  nonperiodic  boundaries. 
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Fill.  13.  The  triple-forward  elimination  ot  the  CCD  Coefficient  matrix  for  EDE  \%iih  nonperiodic  boundaries 
Here  E  denotes  eliminated  coefficients. 


unknowns  at  the  righi  boundary  (x\  .  i ). 


Solving  this  set  of  three  algebraic  equations,  we  obtain  ( 5/  /S.\)\  .  \  and(<52//$.Y:K.:. 

The  substituiion  procedure  starts  from  the  second  right  point  (.v.v  ).  The  dependent  variable 
and  its  first-  and  second-order  differences  at  any  grid  point  (.v, )  tire  computed  from  the 
following  substitution  (/  =  ,\  ,  V  -  1 . 1): 


Appendix  7:  Periodic  CCD  FDK  Solution 

The  structure  of  the  periodic  CCD  PDF.  matrix  is  shown  in  Fig.  14.  We  can  use  a  similar 
triple-forward  elimination  and  triple-backward  substitution  procedures.  Figure  1 5  shows  the 
structure  after  the  triple-forward  elimination  procedure,  where  the  shadowed  areas  mean 
the  eliminated  elements. 


FIG.  14.  Structure  of  tbc  CCD  coefficient  matrix  for  FDF.  with  periodic  boundary 


FIG.  15.  The  triple-forward  elimination  ol  the  CCD  eoefticient  matrix  for  FDF  w  ith  periodic  boundaries. 
Here  □  denotes  eliminated  coefficients 


